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Large-eddy simulation of stratocumulus-topped 
atmospheric boundary layers with dynamic 
subgrid-scale models 

By Inane Senocak 


1. Motivation and objectives 

Earth’s climate and its geographical variation is strongly influenced by cloud coverage. 
It is estimated that about 50% of the earth is covered by clouds at any given time, 
providing a shield from solar radiation. Radiative energy transfer and its interaction 
with clouds play an important role in the thermal structure and stratification of the 
atmosphere. For instance, clouds have high reflectivity in the visible wavelengths, thus 
providing relative cooling of the atmosphere. They also absorb strongly in the infrared 
wavelengths, resulting in heating of the atmosphere (Salby 1996). 

Condensation is the major physical process that is responsible for cloud formation. 
Clouds can be classified into four broad categories, namely: cumulus, cirrus, nimbus and 
stratus (Rogers & Yau 1989). Many other classifications can be derived from combina- 
tions of these four broad categories. A comprehensive description can be found in Scorer 
& Wexler (1967). Among various types of clouds, marine stratocumulus clouds have re- 
ceived increased attention because of the important role they play in the global radiation 
budget. Marine stratocumulus clouds cover about 25% of the Earth’s ocean at any in- 
stant. These are low-level clouds that exist below 1.5 km with several hundred meters in 
thickness and they rarely produce precipitation. Their horizontal coverage is extensive 
and more homogeneous than other type of clouds. Their appearance is grey and a wavy 
undersurface is typical (Kantha & Clayson 2000; Heymsfield 1993; Mason 1975). 

The structure of the marine stratocumulus cloud-topped atmospheric boundary layer 
is driven by both cloud-top radiative cooling and positive buoyancy flux from the surface 
that maintains the atmospheric boundary layer in a well-mixed turbulent state. Above 
the cloud layer, negative buoyancy flux has a stabilizing effect, suppressing the turbulence 
there (Kantha & Clayson 2000). The entrainment of dry air from above the cloud layer 
induces evaporative cooling and entrained air parcel can sink further down enhancing the 
turbulent mixing within the clouds, known as cloud top entrainment instability (Deardorff 
1980a). Clearly, the structure of cloud-topped atmospheric boundary is more complicated 
than a cloud-free atmospheric boundary layer due to strong interactions among cloud 
microphysics, turbulent motions and the radiative energy transfer. 

In large-eddy simulation (LES), three-dimensional, large unsteady flow structures are 
resolved and the effects of the unresolved scales are modeled. A filtering operation is 
applied to the governing equations to distinguish between the resolved scales that are 
computed and smaller scales that are modeled. LES has been widely applied to simulate 
atmospheric boundary layers, partly because of the difficulties involved in observational 
studies and field experiments (Stevens & Lenschow 2001). An extensively studied topic 
is LES of cloud-free convective atmospheric boundary layers. The salient features of 
these boundary layers have been compared to observations and well documented (Kantha 



70 


I. Senocak 


& Clayson 2000). The presence of clouds complicates the problem due to additional 
physics and reliable numerical simulations of cloud-topped boundary layer is still an 
active research area. In the following, several representative studies are briefly mentioned 
to help describe the current state of the knowledge. 

An early study on three-dimensional modeling of cloud-topped atmospheric boundary 
layer is by Deardorff (1980&). He studied the structure of turbulence and entrainment 
within stratocumulus layers with and without cloud-top radiative cooling and suggested 
that simulations need high resolution at the inversion. It was also found that a func- 
tional dependence on Richardson number helps correlate the entrainment rate. Deardorff 
(1980a) has defined a criterion for the cloud top entrainment instability. It was found 
that entrainment rate increases decisively when the equivalent potential temperature gra- 
dient across the cloud top drops below a critical value. It was also shown that a strong 
instability can lead to stratocumulus breakup leading to a scattered cumulus layer. Tag 
& Payne (1987) indicated that in addition to Deardorff’s criterion, the vertical motions 
should exceed a threshold for the cloud breakup to occur. 

Moeng (1986) studied the structure of a stratus-topped boundary layer using LES. It 
was found that the vertical component of the turbulent kinetic energy is generated by 
buoyancy and a portion of this energy is redistributed in the horizontal directions due 
to pressure effects. It was also shown that turbulence is generated more effectively by 
surface heating than cloud-top cooling. 

Bohnert (1993) tested the dynamic procedure for LES of cloud-topped boundary layer. 
Simple parametrization for cloud microphysics and radiation were adopted. The dynamic 
model results were compared to the results of SGS model that were optimized in an ad- 
hoc fashion. The results were found to be comparable. The importance of SGS modeling 
in predicting cloud breakup was also highlighted. 

Stevens et al. (2000) investigated the dependence of an LES model on mesh resolution, 
numerical schemes and SGS model. They provided simulations of varying resolutions for 
the simulation of stratocumulus topped marine boundary layer. Different SGS models 
and advection schemes were tested. It was found that thickness of the inversion layer, 
depth of entraining eddies and shape of the vertical velocity spectra is influenced by 
mesh resolution. Motions at the inversion were found to be underresolved even for the 
finest resolution. The entrainment rate was found to depend on both numerical and SGS 
dissipation. 

Duynkerke, Zhang & Jonker (1995) performed an observational study to describe the 
microphysical and turbulent structure of stratocumulus observed during the Atlantic 
Stratocumulus Transition Experiment (ASTEX). The turbulence kinetic energy budget, 
velocity and temperature variance, and vertical fluxes were calculated. The entrainment 
was found to be very efficient, which resulted in reduction of turbulent kinetic energy 
production due to buoyancy. It was also shown that water vapor flux, liquid water flux, 
and drizzle rate have the same magnitude. 

Stevens et al. (1998) presented an LES study of the ASTEX case. They adopted a 
drop-size resolving cloud microphysics model that enabled them to perform simulations 
with and without drizzle. It was found that inclusion of drizzle in modeling leads to 
sharp decrease in entrainment and turbulent kinetic energy generation by buoyancy. The 
authors have hypothesized that shallow, well-mixed radiatively driven stratocumulus 
cannot persist in the presence of heavy drizzle. 

Duynkerke et al. (1999) did a comparison of actual ASTEX observations with com- 
putations obtained from LES and one-dimensional single column models. The buoyancy 
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flux obtained from LES agrees well with the observations. The authors concluded that 
drizzle has small influence on the buoyancy flux, although significant uncertainty exists 
in its parametrization. 

The objective of the present study is to evaluate the dynamic procedure in LES of 
stratocumulus topped atmospheric boundary layer and assess the relative importance of 
subgrid-scale modeling, cloud microphysics and radiation modeling on the predictions. 
The simulations will also be used to gain insight into the processes leading to cloud top 
entrainment instability and cloud breakup. In this report we document the governing 
equations, numerical schemes and physical models that are employed in the Goddard 
Cumulus Ensemble model (GCEM3D). We also present the subgrid-scale dynamic pro- 
cedures that have been implemented in the GCEM3D code for the purpose of the present 
study. 

2. Numerical formulation 

The numerical model used in this study to simulate cloud-topped atmospheric bound- 
ary layers is the Goddard Cumulus Ensemble model(GCEM3D). Its main features are 
described in Tao & Simpson (1993), Simpson & Tao (1993) and Tao et al. (2003). 

2.1. Governing equations 

Acoustic waves are part of the solution of the compressible Navier-Stokes equations and 
very small time steps are needed to resolve them. On the other hand, acoustic waves do 
not impact the dynamics of thermal convection for low Mach number flows. Therefore, the 
governing equations of motion are simplified by filtering out the sound waves from them. 
The resulting set of equations is known as the anelastic equations (Ogura & Phillips 1962). 
In deriving these equations the basic assumption is to decompose the thermodynamic 
state variables into a horizontally averaged base quantity, which only depends on altitude, 
and a perturbation quantity, which depends on both time and space, as follows 

p(x,y,z,t) =p 0 (z) +p'(x,y,z,t), 
p(x,y,z,t) = p a (z) + p'(x,y,z,t), 

T(x, y, z, t ) = T 0 (z) + T'(x, y, z , t), (2.1) 

where p is the pressure, p is the density of moist air and T is the temperature. The 
following moist equation of state is used in deriving the momentum equations 

p = pRT(l + 0.61q v ), (2.2) 

where R is the gas constant for dry air and q v is the mixing ratio of water vapor. The 
pressure is nondimensionalized according to a reference pressure ( p a ) value 

7i-=( — ) R/Cp , (2.3) 

Po 

where C p is the specific heat of dry air at constant pressure. The potential temperature 
is defined as 

0 =—. (2.4) 

7 r 

The anelastic form of the filtered continuity equation reads as follows 

d(poUj) = 
dxi 


(2.5) 



72 


I. Senocak 


The filtered momentum equations, using the f-plane approximation, are written as 
diii 1 d (pc 


dt 


Po 


~ ~ Cpd °^k + + 9li k + °- 61 ^ _ qi) + Fcoriolis ’ 


( 2 . 6 ) 


where gi is the gravitational acceleration. The Coriolis term that appears in the momen- 
tum equations is written as 


Fcoriolis = [2usin<l>U2, -2wsm^ui,0], 


(2.7) 


where lo is the angular velocity of the Earth and (j> gives the latitude. Tj, represents the 
subgrid scale stress tensor and its particular form is described in section 2.3. The primed 
variables represent deviation from horizontally averaged quantities. 

The equations for potential temperature 6 and the water vapor mixing ratio q v are 
written as follows 


89 1 d(p o 0ui ) 

dt po dxi 



F, 

c z 


(c-e c - e r ) + (f s + f g - m a - m g ) 


C, 


T . , {di ce Sice) T Qradi 

Op 


( 2 . 8 ) 


dg v 1 d(p 0 q v Ui ) 
dt po dxi 


dlT 

dxi 


+ (c - e c 


Cr) T {dice 


(2.9) 


where L v , Lf and L s are the latent heats of condensation, fusion and sublimation, respec- 
tively. The quantities c,e c ,e r ,f,m,di ce ,Si ce represent the rates of condensation, evapo- 
ration of cloud droplets, evaporation of rain drops, freezing of rain drops, melting of snow 
and graupel/hail, deposition of ice particles and sublimation of ice particles, respectively. 
The particular forms for these phase change rates are not explicitly formulated. Instead, 
a saturation adjustment scheme is adopted that calculates the amount of phase change 
rate in order to remove any supersaturated vapor and/or subsaturation of cloud water. 
The saturation adjustment scheme is described in Soong & Ogura (1973) and Tao, Simp- 
son & McCumber (1989). 7 f is the subgrid scale flux of the scalar, which is explained in 
section 2.3. Q ra d is the source due to radiative heat transfer. It is described in section 
2.4. 


2.2. Cloud microphysics model 

The formulation of the cloud microphysical processes is based on solving scalar transport 
equations for each hydrometeor species. The transport equations for cloud water {q c ), 
cloud ice {qi ce ), rain {q r ), snow {q s ), graupel/hail (q g ) are written as 


dg c 1 d{p 0 q c Ui) 

dt po dxi 


dir 

dxi 


+ (c 


e c ) + T qc , 


( 2 . 10 ) 


dq ice 

~dF + 

dg r 1 d{p 0 q r Ui) _ 

dt po dxi 


dg s 1 d{p 0 q s Ui) 

dt p 0 dxi 


1 d{p 0 qiceUi ) d-y f 


Po dxi dx 

1 d{p 0 q r U r ) di^ r 


T {dice Sice) T F, 


Qice "> 


Po dx 3 dxi 

1 d{p 0 q s ll s ) dil 


dx 3 


dx 


+ {m s + m g - f s - f g - e r ) + T qr , 
T {ds s s rn s T f s ) T Tq s , 


( 2 . 11 ) 

(2.12) 

(2.13) 


Po 
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Po 


dx.i 


dxo, 


Ox; 


+ {d g — s g — m g + f g ) + T) 


9s’ 


(2.14) 


where V r , V s and V g are the fall speeds of rain, snow and graupel, respectively. Their 


values are obtained from parameterizations. is the subgrid-scale flux of the scalar. T ^ 
represents the microphysical transfer rates between the hydrometeor species. A total of 27 
different processes are considered and the reader is referred to Lin, Farley & Orville (1983) 
for details of their formulation. To illustrate the parameterizations of these processes, only 
the microphysical transfer rate of rain T qr is described in this report. For instance, if the 
temperature is above 0°C, then the production term for rain is given by the following 
equation 


Tq r = Pracw + Praut + Psacw + Pgacw — Pgmlt — Psmlt + Prevp(1 — Si)- (2.15) 


The terms on the right hand side of the above equation are the accretion of cloud water 
by rain, autoconversion of cloud water to form rain, accretion of cloud water by snow, 
accretion of cloud water by graupel, melting of graupel to form rain, melting of snow to 
form rain, evaporation of rain, respectively. The accretion of rain by cloud water Pracw 
is written as 


TTE rw n 0r aq c T(3 + b) 
PRACW = 


(2.16) 


where E rw is the collection efficiency, which is assumed to be 1, no is the intercept 
parameter of the rain drop size distribution, T is the gamma function, X r is the slope 
parameter in rain size distribution and a, b are empirical constants. Exponential size 
distributions are assumed for rain, snow and graupel/hail. The explicit forms of other 
transformation rates are documented in Lin et al. (1983). 


2.3. Subgrid-scale turbulence models 

Three turbulence closure models will be considered for comparison. These are the subgrid- 
scale kinetic energy model of Klemp & Wilhelmson (1978), which is the base turbulence 
model in the GCEM3D code, the dynamic Smagorinsky model of Germano et al. (1991) 
and the localized dynamic Smagorinsky model of Piomelli & Liu (1995). The imple- 
mentation of dynamic turbulence models by Kirkpatrick (2002) has been coupled to the 
GCEM3D code. In the following sections these models are briefly explained. 


2.3.1. Subgrid-scale kinetic energy model 

A transport equation for subgrid-scale turbulent kinetic energy is solved, which is then 
used to specify the eddy viscosity. The influence of buoyancy on the turbulent motions 
is also modeled. The equation for subgrid-scale turbulent kinetic energy is written as 


dE 1 d{p 0 Eui) 
dt p 0 dxj 


f) * f)ni ■ 

gw\— + 0.61 q’ v - q[) - T ig ^ 


dxj 


(2 ' 17) 


A = (AxAyAz) 1 / 3 , 

K m = C m AE 112 , 

n, - lES, - Km (|| + g). (2.18) 

where K rn is the turbulent eddy viscosity, A is the filter width. The empirical constants 
C e and C m have the values of 0.7 and 0.2, respectively. The subgrid-scale scalar fluxes 
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are modeled as 


= K . ^ 
x ~ h d Xi ' 


(2.19) 


where Kh/K m = 3 is used. The buoyancy flux in a saturated area is computed as 

88 , 


, 8 ' 


'(- + 0.61 g (, - q[) = + I< h 


8q t 

dz 


(2.20) 


where 


A = 


1 1 + 


1.61e Lq v 
RdT 


ft i i eL 2 q v 
° 1 + CpRdT 2 


( 2 . 21 ) 


and e = 0.622. 

In an unsaturated area, the buoyancy flux is computed as follows: 


u>'( j + 0.61 q' v - q[) = -K h ( + 0.61^). 


'8 dz 


dz 


(2.22) 


2.3.2. Dynamic Smagorinsky model 

Smagorinsky model (Smagorinsky 1963) is commonly used in LES to model the subgrid 
scale stresses. It is based on eddy viscosity assumption and can be written as 


nj - \dijTkk = -2CA 2 |S'|5 J : i , 


where C is a dimensionless parameter and IS) = \J 2SijSij . 
defined as 


1 , dm 


dui 


lJ 2 dxd dxi 


(2.23) 

The filtered strain rate is 

(2.24) 


Germano et al. (1991) have introduced the dynamic procedure that offers advantages 
over the original Smagorinsky model. In this model the parameter C is computed at each 
time step from the information already available in the resolved velocity field. The basic 
formalism behind the dynamic Smagorinsky model is explained in the following. 

The dynamic model involves filters of different sizes. In addition to the grid filter, a 
test filter is introduced. The subgrid scale stress tensor based on the grid filter and the 
test filter are written respectively as 


Ti g — UiUj Uj/LL 


Jjl Ujj , 


7 1 j Ui Uj Ui'll j , 


(2.25) 


where the symbols overbar and the hat represent the grid and test filtering operations, 
respectively. Applying the test filter to and subtracting it from 7), yields the following 
identity (Germano, 1992) 


Lij — 7 ’,j Tij — UiUj Ui’llj. 


(2.26) 


The significance of this identity is that it can be computed from the large eddy field. 
Germano et al. (1991) have utilized this identity to dynamically compute the coefficient 
of the Smagorinsky model as follows. 


k J ij ^dijLkk — £%ijC /3ijC, 


(2.27) 



where 
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otij = — 2A 2 SSij, 

fa = -2 A 2 SS ij . (2.28) 


For atmospheric boundary layers, where the horizontal directions are assumed to ho- 
mogenous, the filtering operation is applied only in the horizontal direction and C is 
assumed to be independent of the homogeneous directions and taken out of the filtering 
operator. Equation (2.27) is rearranged to the following form 

Lij — -5ijL kk = CMij, (2.29) 

where 

Mij = ctij - Pij. (2.30) 

Following the method described in Lilly (1992), the coefficient C is computed so as to 
minimize the sum of the squares of the residuals of equation (2. 29). The numerator and 
the denominator are averaged over the horizontal (x,y)-plane. 


( MjjLjj) X y 

( MklMkl)xy 


(2.31) 


Once C is calculated, the subgrid scale stress tensor, given in equation (2.23) is computed. 
In a similar approach, the subgrid-scale flux of a scalar can be computed dynamically 
Moin et al. (1991). If we consider the following eddy diffusivity subgrid-scale model 


vt dcj) 
Pr t dxi ’ 


(2.32) 


where the kinematic eddy viscosity is computed with the dynamic Smagorinsky model 
as v t = (7A 2 |5|. The dynamic procedure can also be applied to compute the turbulent 
Prandtl number Prt . The subgrid-scale scalar flux based on the test filter scale is written 
as 


Pr t OXi 

The test and grid scale fluxes are related to each other by the following identity 

„ „ „ 7 ~ ^A 2 |^| 6>0 A 2 (§f dcj) 

Pi = Gi - 74 = 4 >u.i - <t>Ui = -C — — 1 - C 


Pr t dxi 

The above equation can be recast into the following form 

C 

Pi = -—Ri. 

Pr t 


Prt dxt 


(2.33) 


(2.34) 


(2.35) 


Following the suggestion of Lilly (1992), a least squares procedure is applied to compute 
the value of Pr t . Note that the value of C is determined earlier. 


1 = 1 (PjRi)xy 

Prt C (RrnPm) xy 


(2.36) 


After calculating the Pr t , equation (2.31) is used to compute the subgrid scale flux of a 
scalar. 
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2.3.3. Approximate localized dynamic Smagorinsky model 

The dynamic Smagorinsky model is not general for flows with no homogeneous di- 
rection, because of the need for averaging in the homogeneous directions. Ghosal et al. 
(1995) have addressed the mathematical inconsistencies and proposed a new formulation 
for the dynamic procedure that makes it applicable to arbitrary inhomogeneous flows. 
They have also provided formal justifications for the ad-hoc procedures that have been 
adopted in the early versions of the dynamic model. 

The mathematical inconsistency in previous versions comes from the fact that space 
and time dependent coefficient C is simply taken out of the filtering operation. In the 
variational formulation of Ghosal et al. (1995), C is kept inside and an integral equation 
needs to be solved at each time step to determine C. This method is referred to as the 
dynamic localization model. 

The solution of an integral equation is costly. Piomelli & Liu (1995) have followed a 
simpler approach and proposed the approximate localized dynamic model. In the follow- 
ing this method is briefly described. 

Equation (2.27) is recast in the following form. 

—Cotij = (2.37) 

where 

Lij = Lij — T^dijLkk- (2.38) 

An estimated value C* is assumed, which is the value of C from the previous time step. 
Along with this approximation, a least squares procedure is applied to compute C as 
given below 

(L“- + C*/3ij)cxij 

C = --AA \ (2.39) 

^mn^mn 

The same iterative idea applies to equation (2.31) to compute the Pr t of the subgrid- 
scale flux of a scalar quantity. 


2.4. Radiation model 

In an effort to describe the physical models in GCEM3D code in a single document, we 
briefly summarize the basic formalism of the radiation modeling. 

The plane parallel assumption is typically adopted to model the radiative energy trans- 
fer. The convenience of the assumption comes from the fact that the properties of the 
atmosphere vary sharply with height due to its stratification, hence, the medium is re- 
garded as horizontally homogeneous (Salby 1996). 

The intensity of a radiation pencil, traversing a distance ds along the direction of its 
propagation, changes due to absorbtion, scattering and emission, which can be described 
by the following equation (Liou 2002). 

dl\ = —K\pl\ds + jxpds, (2.40) 

where k\ is the mass extinction cross section due to absorbtion and scattering and j \ 
is the source function due to emission and scattering. In plane-parallel assumption, it is 
convenient to define an optical thickness as 

nO O 

r — j k\dz, dz = fids, fi = cos 9, (2-41) 

J z' 

where z is the vertical direction and 9 is the angle between the path of radiation and the 
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vertical, which is specifically called the zenith angle. The basic equation, describing the 
radiative energy transfer in plane-parallel atmospheres is written as 

M d ^ = j (t; M, 0) ~ J{t\ M, 4>)- (2-42) 

The source function J is defined as 

m', <t>') p (V; </>; m', <t>')dp! d(j)' + 

-MO, ^o)e- T//i ° - (1 - w)B[T(r)], (2.43) 

where P(p, <f>', p' , <j/) is the phase function, which gives the angular distribution of scat- 
tered energy as a function of direction, (f) is the azimuthal angle, Co is the single scat- 
tering albedo, B\T(r)] is the Planck’s function representing the blackbody emission of 
the medium, and F„ is the solar irradiance at the top of atmosphere. Upon solution of 
equation (2.43), the flux density F\ and the total flux density F are computed based on 
the following definitions 

pi pOO 

Fl l (r) = 2 tt / /^(r, T //)//<//, ,F = / F\d\. (2.44) 

Jo Jo 




The heating rate due to radiation that appears in the potential temperature equation, 
Qradi is then computed as 



1 dF(z) 
p 0 C p dz 


(2.45) 


Equation (2.43) is an integrodifferential equation and its numerical solution is quite 
involved due to sharp variation of the atmospheric properties with height. The compu- 
tation of radiation fluxes involves spectral, vertical and directional integrations. Because 
the absorbtion coefficient varies sharply with wave number, the spectral integration is the 
most CPU intensive part. However, the major difficulty comes from the dependence of 
the absorbtion coefficient on pressure and temperature. Hence, effective parametrization 
is an important part in the numerical solution of equation (2.42). Different approaches 
are adopted depending on the nature of the radiation problem. A detailed account of 
atmospheric radiation modeling can be found in Liou (2002). 

The radiative transfer model that is incorporated into the GCEM code is documented 
in Chou & Suarez (1996a, b ). The radiation scheme can model the absorbtion due to water 
vapor, C0 2 , O3, and 0 2 , and scattering by clouds, aerosols and molecules. Fluxes are 
integrated almost over the full spectrum, ranging from 0.175 pm to 10 pm.. The radiation 
field is divided into three distinct regions. In the ultraviolet and photosynthetically active 
regions, the spectrum is divided into 8 bands, in which single O 3 absorption coefficient 
and Rayleigh scattering coefficient is used for each band. The infrared region is divided 
into three bands and the k-distribution method is applied with ten absorption coefficients 
in each band. The <5-Eddington approximation is used to compute the reflection and 
transmission of a cloud and aerosol-laden layer. Two-stream adding method is then used 
to compute the fluxes. 


2.5. Numerical schemes 

In GCEM3D code, the governing equations are discretized on a staggered grid. The mo- 
mentum equations are solved using the second order central difference scheme for the 
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spatial derivatives and the second order accurate leap-frog scheme for temporal deriva- 
tives. A time smoother is adopted to avoid the time splitting problem associated with 
the leap-frog scheme. Forward time differencing and the multi-dimensional positive def- 
inite advection transport algorithm (MPDATA) of Smolarkiewicz & Grabowski (1990) 
are used to solve the scalar transport equations. A direct solver (FFT) is utilized for the 
solution of the pressure Poisson equation. The GCEM3D code has been parallelized to 
run on the SGI Origin clusters using OpenMP (Jin et al. 2003). 

Lateral boundary conditions are periodic and free slip boundary conditions are imposed 
on the top and bottom boundaries for all the variables except the vertical velocity that 
vanishes at the top and bottom boundary. 


3. Future work 

The sugrid-scale turbulence models described in the previous section will be tested 
with progressively refined mesh resolutions, adopting the ASTEX field experiment as 
the test case. The relative influence of subgrid-scale modeling, cloud microphysics and 
radiation modeling on the predictions will be investigated. The simulations will also be 
used to gain insight into the processes leading to the cloud top entrainment instability 
and cloud breakup. Additionally, the dynamic modeling subroutine will be parallelized 
to speed up the overall computation. 

Preliminary computations produced comparable results of wind speeds, soundings of 
cloud water and potential temperature among the subgrid-scale models considered. How- 
ever, the computations also indicated problems with the top boundary resulting in spuri- 
ous cloud formation and unstable stratification. A possible cure to the problem might be 
to use a Rayleigh absorbing layer, which is commonly used in atmospheric simulations, 
in the proximity of the top boundary to damp the vertically propagating gravity waves. 

In this report we have briefly documented the governing equations and physical mod- 
els that are employed in the GCEM3D code. We have also described the subgrid-scale 
dynamic procedures that have been introduced to the GCEM3D code in this study. We 
are in the process of double checking the implementation and the results. The findings 
of the present study will be published in the open literature. 
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